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Introduction 


The pseudospectral differencing methods, involve the exact differentation 
of interpolants which are based on different sets of selected points; each 
method is usually named after the base functions used to expand such 
interpolants. 

We discuss the pseudospectral Fourier and Tchebyshev differencing methods 
— the two most extensively used among all of the above, e.g., the survey of 
Gottlieb, Hussaini and Orszag [5] and the references therein. This stems from 
the possibility of implementing the FFT in these cases: one can efficiently 
travel between the '’physical'’ and 'phase' spaces, making the (global) 
pseudospectral calculations in these two cases, almost as economical as the 
(local) finite difference ones. The definitive advantage of the former lies, 
however, in their remarkable accuracy properties, which is the topic of this 
paper. 

As is well known, the pseudospectral differencing of (sufficiently) smooth 
functions, enjoy 'infinite' order of accuracy; that is, measured w.r.t. the 
inverse number of selected points, the error committed is bounded by any fixed 
polynomial order (e.g., Kreiss and Oliger [8] for the Fourier case, and a 
different detailed study of Canuto and Quarteroni [1], which includes, among 
others, the Tchebyshev case). 

Here we show, that if the function under consideration is further assumed 
to be analytic, then the asymptotic decay rate of the error with either the 
Fourier or Tchebyshev differencing is, in fact, exponential. This should be 
compared with the polynomial decay rate obtained by finite difference/finite 
element differencing methods. 


I 
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We start in Section 2, discussing the Fourier differencing of smooth 
functions: following [8], we first derive the aliasing relation, which 

implies 'infinite' order of accuracy in this case. In Section 3, we show the 
exponential decay rate of the error, with Fourier differencing of analytic 
functions. Tchebyshev differencing method is likewise treated in Section 4: 
after putting the aliasing relation in an identical form to the one obtained 
in the Fourier case, the various error estimates follow along the same lines. 

/ 

Similar to our treatment of the stability question in [15, Part II], we 
emphasize here the central role played by the aliasing relations, from which 
we derive all the results below. Thanks to these aliasing relations, the 
decay behavior is essentially due to the corresponding decay of either 
the Fourier or Tchebyshev coefficients; an exponential decay of the latter is 
widely known in the analytic case. Also, by considering the 

Fourier/Tchebyshev coefficients, the above derivation may still offer an 
exponential decay rate of fractional order in non - analytic, smooth cases 
(e.g., standard cut-off functions). 

In closing, we would like to point out that the above results are 
intimately related to Bernstein s theorem, regarding the exponential 

convergence of best polynomial approximations. Specifically, given an 
analytic function, Bernstein s proof verifies the exponential convergence of 
its truncated Tchebyshev series expansion, e.g., [11, Section 6]. Using the 
Gauss-Tchebyshev rule to compute that expansion's coefficients, we are then 
led to the Tchebyshev interpolant; the further error inferred by such 
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discretization (which is exactly an aliasing error), is known to be also 
exponentially small, e.g., [3, p. 239]. In other words, we conclude that the 
above Tchebyshev interpolent — — so called near minimax polynomial, approxi- 
mates a given analytic function within an exponentially decaying error. In 
fact, the results below indicate that given an analytic function, both the 
Fourier and Tchebyshev interpolants approximate the function and its 
derivatives, within an exponential accuracy. Indeed, these results manifest 
themselves for example, in the global error behavior of pseudospectrally 
solved PDE's, e.g., [5,6,13]. 

Acknowledgement 

I would like to acknowledge Y. Maday for helpful comments concerning this 
work. 


2. Fourier Differencing of Smooth Functions 

Let w(x) be a 2x-periodic function, whose values, w^ = w ( x v )» are 
assumed known at the 2N equidistant grid points x v =vh, h=^, v=0,l,*»*2N-l. 
The (pseudospectral) Fourier differencing of such function, refers to 
differentation of the trigonometric interpolant of these grid values: one 

constructs the trigonometric interpolant^^ 

N 2N-1 

(2.1) w(x) = w(x;N) = l w e ipx , w p = -tf* l w v e" ipvh , 

p=-N v=0 


(1) (double) primed summation indicates halving first (and last) terms* 
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and use its derivative 


dw/ \ _ * px V 

^(x v ) - l ipw p e 


to approximate the 'true' value, - r -{x=x ). 

dx v 

In order to examine the error we commit by such approximation, it is 
convenient to work with Sobolev space W s , defined for integral orders s, 


( 2 . 2 ) 


w s = W® = {w(x)| Jwj 2 = l j — k — | 2 < -»}, 

W k=0 dx L [0 ,2tt] 


and extended by interpolation for fractional orders. Thanks to Plancherel's 
formula, W s is isometrically isomorphic to H s : assuming w(x) admits a 
formal Fourier expansion 


(2.3a) w(x) ~ £ w(p)e ipx , w(p) = w(?)e“ ip ^d5; 


then we can equally work with H s , s real, which consists of those functions, 
w(x), having a finite Sobolev norm of order s, 


(2.3b) H s = {w(x)| jw| 2 s = l (1+| p | ) 2s | w(p) | 2 < «} 


The following lemma, relating the Fourier coefficients of w(x). 


(2.4) 


Kp) = / w(5)e“ ip ^d?, 

5=0 


< P < 
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with those of its trigonometric interpolant, w(x;N), 

2N-1 

(2.5) w p = I w(x v )e“ i P vh , -N <_ p <_ N, 

v=0 

is in the heart of our discussion (e.g. , Kreiss and Oliger [8]). 

Lemma 2.1. (Aliasing) . 

Assume w(x) is in H s , s > V 2 • ^ 2 ^ Then the following equality holds 

00 

(2.6) w p = I w(p+2kN), -N < p < N. 

k=-» 

Proof : It is well known that the Fourier expansion in (2.3a) converges in 

this case, e.g., [16, Chapter II] ; inserting that expansion, evaluated at 
x = x^, into (2.5), we obtain 

% 

2N-1 2N-1 » . 

(2.7) ip --Jr I w(x v ).- 1 " Uh - 1- I [I i(q)e q v ]e- 1 I>' ,h . 

v=0 v=0 q=-« 

By Cauchy-Schwartz inequality the inner summation is absolutely convergent 

(2 ‘ 8) l |w(q)| < [ I 7T~T7?] 7 *|w| , 

q=-oo q=l ( ) “ 

hence summations on the right of (2.7) can be interchanged. By so doing, the 
desired result follows 


( 2 ) 


This smoothness assumption on 


w(x) 


can be relaxed. 


T 
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» 2N-1 « 

w p » l w(q)*^* l e iv ^ q “ p ^ h = l w(p+2kN) , 
q=-°° v=0 k=-°° 

noting that the second summation in the middle term vanishes unless q — p = 
0(mod 2N), i.e., q ■ p + 2kN. 


Equipped with the aliasing lemma, we now may turn to estimate the error 
between w(x) and its equidistant interpolant w(x) : rewriting 

w(x) = [ l + l ]w(p)e ipx , 

|p|<N |p|>N 

and, with the help of (2.6), 


w(x) = l w(p)e lpx + l [ l w(p+2kN)]e ipx , 

| p | <N |p|<N k#0 

the difference w(x) - w(x) is readily verified to equal 

(2.9) w(x) - w(x) = - l [I w(p+2kN)]e ipx + l w(p)e lpx . 

|p|<N k^O |p|>N 

The first summation on the right represents aliasing of the higher modes with 
the lower ones, |p| N, while the second summation consists of the truncated 
higher mode, |p| >_ N. A quantitative study of both terms gives us (compare 
e.g. , Kreiss and Oliger [9], Pasciak [12]) 

Lemma 2.2. (Error Estimate) . 

g 1 

Assume w(x) is in H , s > -j. Then for any real a, 0 < a s, we have 


Proof 


Starting with (2.9), then by definition 


( 2 . 11 ) 


jw(x) - w(x;h)| 2 a = l (l+|p|) 2a | l w(p+2kN)| 2 

H |p|<N k*0 

+ I (l+|p| ) 2CT |w(p) | 2 . 
|p|>N 


Cauchy-Schwartz inequality implies 


| I w(p+2kN) | 2 <_ l (l+|p+2kN| ) 2s * |w(p+2kN) | 2 » £ (l+|p+2kN| )“ 2s , 
k*0 k^O k*0 


with the second summation not exceedng a value of 


l (l+|p+2kN|r 2s < 2N~ 2s * l (2k-l)” 2s , |p | < N. 
k*0 k=l 


Inserted into (2.11), we find that the aliasing part of the error given in the 
first term on the right, is bounded by 


2N _2s * l (2k-l)” 2s * l N 2a l (l+|p+2kN|) 2s |w(p+2kN)| 2 
k=l |p|<N k*0 


< 


CD 

2* l (2k-l)” 2s * 
k=l 


A2(s-o) 

'■N'' 



The truncation error, given in the second term on the right of (2.11), is 
equally found to be bounded by 
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r" N 2(0 "s). (1 +| p | ) 2s |w(p) | 2 <. (fc 2(s ” a) ‘|w|2 . 

|p|>N 


Added together, the last two estimates yield (2.10). 


Remark 1 . Observe that requiring w(x) to have more than "one-half" 
bounded derivative enable us to control the aliasing part of the error; apart 
from that restriction, there is an error in decay in any Sobolev norm weaker 
than that of w(x) , which is equally due to aliasing and truncation errors. 


Remark 2. The aliasing relation (2.5) for the zeroth mode p = 0, implies 
that the trapezoidal rule is highly accurate for the integration of smooth 
periodic functions (Davis and Rabinowitz [3]): indeed, the error committed in 
this case is solely due to aliasing 


2N ~ 2it 

Tf I w(x v ) - Tf-*/ w(5)d£ = I w(2kN) . 
v=0 0 k*0 


This allows us to replace the H a -norm, measuring the error on the left of 
(2.10), with its more applicable discrete counterpart (Gottlieb et.al. [5]): 


|w(x)-w(x;N| | | 


H 


a 2N„ (k) 

l ir l t W - k 

k=0 ZW v=0 dx k V dx* 


d (k) w. m .i 2 
~(x v ;N)J 


a integral. 


Returning to our original question, we find — choosing o = 1 in Lemma 
2.1 — that the error in Fourier differencing does not exceed 


( 2 . 12 ) 


1 < Consf|w| •(^) s -* , 

H 


for arbitrary real s, s > 1. The norm on the left refers, of course, to 


the H® = norm of the error, with a uniform Constant = 2 on the right. It 
can be replaced, in fact, by any other reasonable (possibly discrete) norm: 
for example, Sobolev's inequality implies for the somewhat more applicative 
maximum norm 


Max 

0<V<2N-1 


dw, x dw, 
j — (x ) ~ j — (x ;N) 
dx' v' dx v’ 


< Const 


r*I w l 


H" 


V » 


s > o 


>! 


Consider now a sufficiently smooth 2ir-periodic function w(x). 

Differencing such function by local methods, such as finite difference or 
finite element methods, leads to an error bound of the type (2*12) with a 
finite, fixed^ ^ degree, polynomial decay; the latter is usually identified 
with the accuracy order of the differencing method. With this terminology in 
mind, the (global) Fourier differencing method was thus shown to be 
"infinitely” order accurate: the discretization error decays faster then any 

fixed degree polynomial rate, e.g., [1-2], [4-7], [14-15]. It is worth 

emphasizing that phrasing the error estimate (2.12) as "infinite" order of 
accuracy, is limited on both accounts: 

1. Consider a sufficiently smooth function w(x) in H s , s » 1. The 
error^s order of magnitude for a given Fourier differencing of such 
functions, may be difficult to calculate: an a' priori knowledge 

regarding the size of the factors | w| k s, is required in this case. 


(3) 


That is, independent of N 
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OP 

2. Assume w(x) is a C - function. One cannot detect the exact 

asymptotic decay rate , according to the error estimate (2.12): because of 

its factor dependence on the power s — when s increases so does 

Jwl — one may not conclude, for example, an exponential convergence 
H 

rate simply by placing arbitrarily large powers s, since the optimal s 
depends of course (usually in unknown manner) on N. 


3. Fourier Differencing of Analytic Functions 

In this section, we show that the Fourier differencing of 27r-periodic 
analytic functions, admits an exponentially decaying error; furthermore, in 
some cases, the error's order of magnitude may be calculated as well. 

To this end, assume 

(3.1a) -Hq < Imz < rig 

to be the strip of analyticity where w(z) admits the absolutely convergent 
expansion 

CO 

(3.1b) w(z) = l w(p)e i P z , | Imz | <. n < J1 q. 

p=-» 

Denoting 

(3.2) M(n) = Max |w(z)|, 

| Imz \<n 

we may now state 
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Theorem 3.1. 

Assume w(x) _is_ 2ir -periodic analytic, with analyticity strip of width 
2 Hq. Then for any n , 0 < n < Pq, we have 

I 

(3.3) |g(x) - ^(x;N)| < 4M(p)( Ctg ffi T ^ ) 2 »Ne~ Nn . 

e —1 

Proof . Making the change of variables, ? = e* Z , then v(?) = w(z=-iLog?) 
admits the power series expansion 

00 

(3.4) v(?) = w(z=-iLog?) = j> w(p)? p . 

psr— OO 

By the periodic analyticity of w(z) in the strip |lmz|<riQ, v(t) is found 

-n 0 fin 

to be single valued analytic in the corresponding annulus e <£<e , whose 

Laurent^s expansion is given in (3*4): 

(3.5) w(q) = -jjfx* / d? > e 0 < r < e °. 

kl=r 5 

To estimate the error of Fourier differencing in this case, we employ (2.11) 
with 0=1, obtaining 


(3.6) Jw(x)-w(x;N)j <. N 2 *£ | £ w(p+2kN)| 2 + £ (l+|p| ) 2 |w(p) | 2 

H |p|<N k*0 |p|>N 

Using (3*5), we sum the aliased amplitudes 


[ I + I ] w(p+2kN) - Tfr [ / p+l^^) + , { -i 5 P+I(5-^I1)1» 

k<0 k>0 |c|=r 5 L) |?|=r 1 ^ l) 


i n >l, 


T 
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so that the first term on the right of (3.6) does not exceed a value of 

<3. 7a) 4N 2 "“‘^ r|) l"e 2np < 4M 2 (n) ct8 M Nn) N 2 e~ 2Sl1 . 

(e |p|<N e "- 1 

The truncation contribution to the error in the second term on the right of 
(3.6), does not exceed^^ 

(3.7b) 4M 2 (n)[ f (l+p 2 ).- 2t " > + l (l+p 2 ). 2 ”’ 1 ’] < 8 4^ hV 21 " 1 . 

p>N p<-N e n -l 

Adding the last two bounds, yield (3.3). 

Remark 3 . According to the above derivation, the asymptotic dependence of 
the overall error on N, is due to equal size contributions of both the 
aliasing and truncation parts. However, one can do better with regard to the 
truncation error: indeed, let us denote 

k 

R (n) - e 2kn - I Max |v (J) (5)|j 
11 j-0 |c|-a n 

then by invoking the relation 

qw(q) = / ^k?)C q dC, 

U|=r 

^^We assume N is sufficiently large, N > (e^-1) 
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the truncation contribution in (3.2b) is, in fact, found to be bounded by 


M?(n) 


2M?(n)[ I'e“ 2np + l e 2t1p ] < 4 -111 e“ 2Nri . 
1 p>N p<-N £Tl ’ 


e 2n -l 


Compared with the truncation estimate in (3.7b), we see that the loss of the 
2 

N -factor is regained here. 


Remark 4 . Estimate (3.3) shows that the error with Fourier differencing 
of an analytic function w(x), decays exponentially w.r.t. its asymptotic 
dependence on N; furthermore, equipped with a bound on w(x) when moved into 
the complex plane, one can estimte the size of the error in this case, using 
the somewhat more aesthetic upper bound 


(3.8) 



4M(n) 

— sinh(n) 


■Ne 


-Nn 


Remark 5. The exponential convergence follows for derivatives higher than 
one: with the usual loss of a factor of N for each derivative, we obtain 


O. 9 ) ! w(x)-w(x;N)| < Const ff .- Bl jffi^ N a e' Nn * 

The preferable discrete estimates, follow along the lines of an earlier 
remark, or alternatively, using Sobolev inequality to implement L°° error 


estimates 
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4. Tchebyshev Differencing - The Non-Periodic Case 

In the non-periodic case, the Tchebyshev differencing is usually 
advocated, e.g., [1-2], [4-6], [10], [13-14], Let w(x) be defined for 
”1 <. x < 1, and assume its values w v = w(x v ) are known at the N+l 
gridpoints Xy=cos(vh), h = -^ , v = 0,1,»**N. The (pseudospectral) Tchebyshev 
differencing of such function refers to differentiations of the polynomial 
interpolant of these gridvalues: one constructs the polynomial interpolant 

(4.1) w T (x) = w T (x,N) = 'wT(x), w = i *r'w v Tp<x v ), 

p=0 v=0 

in terms of Tchebyshev polynomials Tp(x) = cos [p(cos”*x) ] , and use its 
derivative 

dw, v dT p < - x 

Hx^ x v) = 4 w ^(x=x v ) 

p=0 

to approximate the 'true' value, -^x=x ). The latter summation, can be 
translated into standard cosine FFT-like summation using a single two-step 
recurssion formula, see [4-6]; thus Tchebyshev differencing admits a fast 
efficient implementation. 

To measure the error in this case, one usually employs the appropriately 
weighted Tchebyshev norm 



} w 2 (x) 

-1 a -* 2 ) 1 ' 2 


dx 


and the corresponding weighted spaces under the W® norm, s integral. 
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(4.2) 


W® = { w(x) | 



Tchebyshev spaces W^ of factional order s are suitably interpreted by 
interpolation. 

We have found it more convenient, however, to work below within the spaces 

g 

Up s real: assuming w(x) admits a formal Tchebyshev expansion 


(4.3) 


V w(p)T (x), w(p) = fj 


1 w(5)T n (5) 


d-r) j 


then, in complete analogy with (2#3b), we introduce 


CO 

(4.4) H| = {w(x)| |w| 2 s = l (l+p) 2s |w(p)| 2 < •}. 

% p=0 

Unlike the Fourier case (endowed with the usual Euclidean 

s s 

weighting), and H T are not equivalent unless s = 0, in which case they 

are in fact isometrically isomorphic by the Tchebyshev transform 


(4.5) 




Making use of the inverse inequalities of Canuto and Quarteroni [1], will 
enable us, later on, to recover the H®-estimates derived below, within the 

g 

more standard W^-spaces. We begin with the aliasing relation in this case, 
which reads (e.g. Gottlieb [4], Reyna [14]) 
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Lemma 4.1. (Aliasing). 


- s X 

Assume w(x) is in IL,, s > ■=■. Then the following equality holds 


(4.6) 


w p - w(p) + £ (w(-p+2kN) + w(p+2kN)], 0 p N. 

k=l 


Proof . Inserting the Tchebyshev expansion in (4.3) evaluated at x = x y , 
into the discrete Tchebyshev coefficient in (4.1), we find 


N «° 

= 2 v~r V 


w p =& l [ I w(q)T q (x v )]T (x v ) = |fw(q)[ 1"t (x v )T (x v )]. 


v=0 q=0 


q=0 v=0 


To calculate the inner summation on the right, we employ the identity 
2 T q (x)T q (x) = T p+q (x) + T | p - q | (x) > while noting that I T j (x v ) vanishes, 
unless j = 0(mod 2N) in which case it equals one. Hence, we end up with 


W p I v/(q)[S + ^qO*^pO + E ^q,2kN±p^ 

q=0 k=l 


and (4.6) follows. 


Let us define T_^(x) = T^(x) so that w(-p) = w(p). Tchebyshev 
expansion (4.3) takes now the Fourier-like symmetric form 


1 _ A 

w(x) ~ l w(p)T (x) 


with an aliasing formula, Identical to the one we had before in Lemma 2.1: 


(4.8) 


w = l w(p+2kN). 
k=— 

Hence, we can equally conclude the corresponding error estimate, which we 
quote from Lemma 2.2. 

Lemma 4.2. (Error Estimate) 

s 1 

Assume w(x) is in H^, s > j. Then for any real a, 0 a < s, we have 

1 

OO o- 

(4.9) | w(x) - w T (x;N)j g < 2(l+2« l (2k-l) 2S ) 2 *|w| S *(^) S ° • 

H t k=l H t 

Setting a = 0 in (4.9) gives us, in view of (4.5) 

1 

(4.10) |w(x) - w t (x;N)| q [2ir(l+2- l (2k-l)“ S )] -|w| 

W T k= 1 H,j, 

Using the inverse inequality [1, Lemma 2.1], one can 'raise' the Sobolev norm 
on the left of (4.10), obtaining (for details see Canuto and Quarteroni [1, 
Theorem 3.1], Maday and Quarteroni [10]) 

Corollary 4.3. (Error Estimate) . 

s 1 

Assume w(x) is in W^, s > • Then for any real a, 0 a s, we have 

~ 1 s “2° 

(4.11) |w(x) - w T (x;N)j g < Const g «|w| s *%) 

w T w T 

Thus, each derivative infers a lost of factor in this case, rather than 

the usual factor N associated with the Fourier differencing. 
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Remark 6. According to Y. Maday (private communication), the factor 

dependence on the right of (4,11) is factorial. Const ~ s!. 

s 

We turn now to consider the case where w(x) is analytic in the interval 
[-1,1], To this end, we employ Bernstein's regularity ellipse, E r , with 
foci ±1 and with sum of its semiaxis equals r, e.g., [11, Section 6], 

Denoting 

(4.12) M_(r) = Max |w(z)|, 

zeE 

we may now state r 

Theorem 4.4 

Assume w(x) is analytic in [-1,1], having a regularity ellipse whose 
sum of its semiaxis equals r Q = e^ > 1. Then for any n, 0 < n < n Q , we 
have 

( 4 * 13 ) |w(x) -w T (x;N)J l < 8M(n)(-- tg ^ Nn) ) 2 »Ne~ Nn . 

e -1 

c+c -i 

Pr °o f » The transformation, -^2 — = z , takes the regularity ellipse E 

-1 r ° 
in the z-plane, into the annulus r Q <|c|<r Q in the 5 -plane. Hence, 

t +c" 

v(?)=2w(z= — ^ — ) admits the power series expansion 

(4.14) v(?) = 2w( z= ^ ■ ~ 1 ] = l w(p)?P, tq 1 < |c| < r 0 =e n ° . 

P— * 

Indeed, upon setting C = e and recalling that w(— p) =* w(p), the above 
expansion clearly describes the real interval [-1,1], 
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w(z=COS0)= £ w(p)cos(p6). 
p=0 

For the Laurent's expansion given in (4.14), we then find 


(4.15) 


w(q) = twt* . { 


“^0 


I c qH 
?|=r « 


< r < e 


n 0 


Comparing (4.15) and (3.5), we end up with the same Cauchy integral formulae 
for the amplitudes in both the Fourier and Tchebyshev expansions; coupled with 
the identical aliasing relations, (4.13) follows along the lines of Theorem 

3.1. 


Remark 7. As before, the factor ( ct gh(N , n) 

e - 1 1 

be replaced by the more aesthetic bound of s j -n h(i^) 


on the right of (4.13), 
-, yielding 


can 


< 4 - 16 > l"« - > < 8 5lfRS) ,,rN "- 

H T 

Next, an exponential error estimate in terms of the Sobolev norm W^, can 
be derived: with the loss of an additional factor of N in the spirit of an 

earlier remark, we then find 


Corollary 4.5. 

Assume w(x) is analytic in [—1,1]. Then we have 


(4.17) |S£(x) 


dw_ 

tj^-(x;N)I t < Const 


M(tQ 

sinh(n) 


2 -Nn 
N e 


y 


0 < n < 


r 
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Higher derivatives can be estimated in a similar manner; in particular, since 

oo 1 /2 + G 1 

L CH C.W , a discrete maximum estimate follows. 


Corollary 4.6. 

Assume w(x) is analytic in [-1,1]. Then we have 


(4.18) 


Max 

0<v<2N- 


_ 1 lS (x v ) " S ( V N) I < Const< 7 i ~ nh(n) N4e ~ Nn * 0 < n < 11 


O’ 


The fourth power of N on the right of (4.17) can be improved to be 
n 5/2 + e > ^ imbedding W* ,# ° in + e . 
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